#include "PairedReadsMapping.h"
#ifdef _OPENMP
#include <omp.h>
#endif
// Given the list of paired read sets in the single file
// format and the index table, this function maps reads int parallel.

int parallelMapPairedReads(vector<string>& readSetsList, GIndexTQ& table1, GIndexTQ& table2, MappingOpts& P)
{
    P.bDiscardReadWithN = false;
    P.clearOutputFileName(readSetsList.size() > 1);
    int i;
#ifdef _OPENMP
    int numberOfCPUs = omp_get_num_procs();
    LOG_INFO("\nInfo %d: %d CPUs detected. %s.\n",\
             INFO_LOG, numberOfCPUs, BLANK_LINE);
    #pragma omp parallel for
#endif
//__OPENMP_FOR_PARALLEL__(#pragma)#ifdef _OPENMP
    for (i = 0; i < (int)readSetsList.size(); i++) {
        CPairedReadsMapping mapping(P);
        const char* readSetName = (readSetsList.at(i)).c_str();
        if (checkFileExist(readSetName)) {
            // assume the paired end reads is in the format
            // that concatenate paired read in 5'-3' and 3' to 5'.
            bool in5to3cat3to5Format = true;
            int allowedNumOfNinRead = P.readLength; // TODO: Don't throw reads with N in the paired-end mapping
            // TODO: Becareful that P.readLength + P.readLength2 may not be set yet
            CPairedReadsSet pairedReadSet\
            (readSetName, P.readsFileFormat, P.readLength + P.readLength2, in5to3cat3to5Format, allowedNumOfNinRead);
            if (P.bIgnoreQS) {
                pairedReadSet.ignoreQScores();
            }
            CReadInBitsSet& readSet1 = *pairedReadSet.F_Reads;
            CReadInBitsSet& readSet2 = *pairedReadSet.R_Reads;
            TIME_INFO(mapping.mapPairedReadsInPairedFiles(readSet1, readSet2, table1, table2), "Mapping takes");
        }
    }
    return(0);
}

// Given two mapped paired read set list and the index table,
// this function maps reads in parallel
int parallelMapPairedReads(vector<string>& readSetsList1, vector<string>& readSetsList2, GIndexTQ& table1, GIndexTQ& table2, MappingOpts& P)
{
    P.bDiscardReadWithN = false;
    ASSERT_TRUE(((int)(readSetsList1.size()) == (int)(readSetsList2.size())), "Read sets are not in paired");
    P.clearOutputFileName(readSetsList1.size() > 1);
    int i;
#ifdef _OPENMP
    int numberOfCPUs = omp_get_num_procs();
    LOG_INFO("\nInfo %d: %d CPUs detected. %s.\n",\
             INFO_LOG, numberOfCPUs, BLANK_LINE);
    #pragma omp parallel for
#endif
//    __OPENMP_FOR_PARALLEL__(#, pragma)
    for (i = 0; i < min((int)readSetsList1.size(), (int)readSetsList2.size()); i++) {
        CPairedReadsMapping mapping(P);
        const char* readSetName1 = (readSetsList1.at(i)).c_str();
        const char* readSetName2 = (readSetsList2.at(i)).c_str();
        if (checkFileExist(readSetName1) && checkFileExist(readSetName2)) {
            if(isLongPairedReads(P.readLength, P.readLength2)) {
                CLongReadsSet rSet1, rSet2;
                initPairedReadSet(rSet1, rSet2, readSetName1, readSetName2, P);
                TIME_INFO(mapping.mapPairedReadsInPairedFiles(rSet1, rSet2, table1, table2), "Mapping takes");
            } else if(isLongShortPairedReads(P.readLength, P.readLength2)) {
                CLongReadsSet rSet1;
                CReadInBitsSet rSet2;
                initPairedReadSet(rSet1, rSet2, readSetName1, readSetName2, P);
                TIME_INFO(mapping.mapPairedReadsInPairedFiles(rSet1, rSet2, table1, table2), "Mapping takes");
            } else if(isShortPairedReads(P.readLength, P.readLength2)) {
                CReadInBitsSet rSet1, rSet2;
                initPairedReadSet(rSet1, rSet2, readSetName1, readSetName2, P);
                TIME_INFO(mapping.mapPairedReadsInPairedFiles(rSet1, rSet2, table1, table2), "Mapping takes");
            }
        }
    }
    return(0);
}

CPairedMappingFlags::CPairedMappingFlags(void)
{
    this->bNoMapping = true;
    this->sepMore = false;
    this->sepLess = false;
    this->pairedOnExpStrand = false;
    // should be constant
    this->bPrintAllMapping = false;
    this->mapSOLiDRead = false;
}

CBestPairedMapping::CBestPairedMapping(void)
{
    this->validMappingNo = 0;
    this->bestMappingNo = 0;
    this->minDiff = MAX_READ_LENGTH;
}

CBestPairedMapping::~CBestPairedMapping(void)
{
}

inline void CBestPairedMapping::update\
(CMappingResult &m1, CMappingResult &m2, bool excludeAmbigousRead)
{
    unsigned int diff = m1.uiDiff + m2.uiDiff;
    if (diff < this->minDiff) {
        this->bm1 = m1;
        this->bm2 = m2;
        minDiff = m1.uiDiff + m2.uiDiff;
        this->bestMappingNo = 1;
    } else if (diff == minDiff) {
        this->bestMappingNo++;
    }
    this->validMappingNo++;
}


CPairedReadsMapping::CPairedReadsMapping(void)
{
    this->initialization();
}

CPairedReadsMapping::~CPairedReadsMapping(void)
{
}

CPairedReadsMapping::CPairedReadsMapping(const MappingOpts P): CReadsMapping( P )
{
    this->initialization();
}

void CPairedReadsMapping::initialization(void)
{
    this->noOfPairsInRange = 0;
    this->noOfSingle1stEndMapped = 0;
    this->noOfSingle2ndEndMapped = 0;
    this->noOfAmbiguousPairs = 0;
    this->noOfPairsSepLess = 0;
    this->noOfPairsSepMore = 0;
    this->noOfPairsSepMoreAndLess = 0;
    this->noOfExpMappedPairedStrand = 0;
}
/* Deprecated
int CPairedReadsMapping::mapPairedReadsInASingleFile\
(CPairedReadsSet& readSet, GIndexTQ& table)
{
    CReadInBitsSet& readSet1 = *readSet.F_Reads;
    CReadInBitsSet& readSet2 = *readSet.R_Reads;
    cout << "Start mapping " << readSet1.InputFile << " and " << readSet2.InputFile << endl;
    getReadsFileFormat(readSet.InputFile, opt.readsFileFormat);
    if (this->setUpIO4Aligment(readSet.InputFile, table) != 0) {
        LOG_INFO("\nInfo %d: Fail to setup I/O files.", ERROR_LOG);
        return(1);
    }
    while (readSet.get_next_capacity_reads_pairs_from_single_file()) {
        mapPairedReads(readSet1, readSet2, table);
    }
    this->printMappedPairStats(cout, readSet1, table.uiSubDiffThreshold * 2);
    this->tearDownIO4Aligment();
    return(0);
}
*/

int CPairedReadsMapping::mapPairedReads(CReadInBitsSet& readSet1, CReadInBitsSet& readSet2,
                                        const GIndexTQ& table1, const GIndexTQ& table2)
{
    CAlignmentsQ& q1 = this->alignmentsQ[0];
    CAlignmentsQ& q2 = this->alignmentsQ[1];

    int bufferedReadNo = checkPairedReadSetSize(readSet1, readSet2);
    vector<CReadInBits>::iterator it1, it2;
    it1 = readSet1.pReadsSet->begin();
    it2 = readSet2.pReadsSet->begin();
    for	(int i = 0; i < bufferedReadNo; i++, it1++, it2++) {
        q1.read = *it1;
        q2.read = *it2;
        readSet1.get_read_id(i, q1.tag);
        readSet2.get_read_id(i, q2.tag);
        q1.qualityScores = readSet1.getQScoresPtr(i);
        q2.qualityScores = readSet2.getQScoresPtr(i);

        this->queryARead(*it1, table1, q1);
        this->queryARead(*it2, table2, q2);
        // statistics and output
        if (q1.load > 0 && q2.load > 0) {
            // TODO check the funciotn
            this->dealMappedPairedReads(table1, table2);
        } else if (this->opt.bPrintUnMappedReads) {
            // TODO to be done
            // dealMissedPairedRead(m1, m2);
        }
    }
    return(0);
}

/* TODO List: for mapping long/short paired read.
 * (1) Test if single SOLiD long read works
 * (2) Test if SOLiD read will work for 2X, X.
 * (3) Make long/short paired read, X>64, Y<64 work.
 * (4) Make long/short paired read, X<64, Y<64 work.
 */
int CPairedReadsMapping::mapPairedReads(CLongReadsSet& readSet1, CReadInBitsSet& readSet2, const GIndexTQ& table1, const GIndexTQ& table2)
{
    /*
    if(isMappingSolidPaired(readSet1, readSet2, table1, table2)) {
        LOG_INFO("\nInfo %d: PerM currently doesn't support SOLiD reads paired with different lengths.", ERROR_LOG);
        return(1);
    }*/

    CAlignmentsQ& q1 = this->alignmentsQ[0];
    CAlignmentsQ& q2 = this->alignmentsQ[1];
    CReadInBitsSet& readSetA1stHalf = *(readSet1.F_Reads);
    CReadInBitsSet& readSetA2ndHalf = *(readSet1.R_Reads);

    int bufferedReadNo = checkPairedReadSetSize(readSet1, readSet2);
    vector<CReadInBits>::iterator itA1, itA2, itB1;
    itA1 = readSetA1stHalf.pReadsSet->begin();
    itA2 = readSetA2ndHalf.pReadsSet->begin();
    itB1 = readSet2.pReadsSet->begin();
    for (int i = 0; i < bufferedReadNo; i++, itA1++, itA2++, itB1++) {
        this->printCheckPointInfo(i);
        CMappingResult m1, m2;
        this->getLongReadInfo(readSet1, i, *itA1, *itA2, m1);
        this->getReadInfo(readSet2, i, *itB1, m2);

        this->queryALongRead(*itA1, *itA2, table1, q1);
        this->queryARead(*itB1, table2, q2);

        // statistics and output
        if (q1.load > 0 && q2.load > 0) {
            // TODO fix this to be LongShort
            this->dealMappedLongPairedRead(table1, table2, m1, m2);
        } else if (this->opt.bPrintUnMappedReads) {
            // TODO to be done
            // dealMissedPairedRead(m1, m2);
        }
    }
    this->iReadCounter += bufferedReadNo;
    return(0);
}

int CPairedReadsMapping::mapPairedReads(CLongReadsSet& readSet1, CLongReadsSet& readSet2, const GIndexTQ& table1, const GIndexTQ& table2)
{
    CAlignmentsQ& q1 = this->alignmentsQ[0];
    CAlignmentsQ& q2 = this->alignmentsQ[1];
    CReadInBitsSet& readSetA1stHalf = *(readSet1.F_Reads);
    CReadInBitsSet& readSetA2ndHalf = *(readSet1.R_Reads);
    CReadInBitsSet& readSetB1stHalf = *(readSet2.F_Reads);
    CReadInBitsSet& readSetB2ndHalf = *(readSet2.R_Reads);

    int bufferedReadNo = checkPairedReadSetSize(readSet1, readSet2);
    vector<CReadInBits>::iterator itA1, itA2, itB1, itB2;
    itA1 = readSetA1stHalf.pReadsSet->begin();
    itA2 = readSetA2ndHalf.pReadsSet->begin();
    itB1 = readSetB1stHalf.pReadsSet->begin();
    itB2 = readSetB2ndHalf.pReadsSet->begin();
    for (int i = 0; i < bufferedReadNo; i++, itA1++, itA2++, itB1++, itB2++) {
        this->printCheckPointInfo(i);

        CMappingResult m1, m2;
        this->getLongReadInfo(readSet1, i, *itA1, *itA2, m1);
        this->getLongReadInfo(readSet2, i, *itB1, *itB2, m2);
        this->queryALongRead(*itA1, *itA2, table1, q1);
        this->queryALongRead(*itB1, *itB2, table2, q2);

        // statistics and output
        if (q1.load > 0 && q2.load > 0) {
            this->dealMappedLongPairedRead(table1, table2, m1, m2);
        } else if (this->opt.bPrintUnMappedReads) {
            // TODO to be done
            // dealMissedPairedRead(m1, m2);
        }
    }
    this->iReadCounter += bufferedReadNo;
    return(0);
}

inline bool isExpPairedMappedStrand(bool firstEndFirst, char end1Strand, char end2Strand, bool bSOLiD)
{
    if (bSOLiD) {
        if (firstEndFirst && (end1Strand == '+') && (end2Strand == '+')) {
            return(true);
        } else if (!firstEndFirst && (end1Strand == '-') && (end2Strand == '-')) {
            return(true);
        }
    } else {
        if (firstEndFirst && (end1Strand == '+') && (end2Strand == '-')) {
            return(true);
        } else if (!firstEndFirst && (end1Strand == '-') && (end2Strand == '+')) {
            return(true);
        }
    }
    return(false);
}

inline bool isValidPaired(CMappingResult& m1, CMappingResult& m2, MappingOpts& opts)
{
    bool allowSameStrand = !(opts.frOnly);
    bool allowDiffStrand = !(opts.ffOnly);
    if (strcmp(m1.RNAME, m2.RNAME) == 0) {
        if(m1.strand == m2.strand) {
            return(allowSameStrand);
        } else {
            return(allowDiffStrand);
        }
    }
    return(false);
}

// TODO becareful about the big separation that overflow to negative value
inline int getSep(unsigned int uiPos1, unsigned int uiPos2, bool expM2gtM1)
{
    int range = 0;
    if (uiPos1 < uiPos2) {
        range = (int)(uiPos2 - uiPos1);
        if(!expM2gtM1) {
            range *= -1;
        }
    } else {
        range = (int)(uiPos1 - uiPos2);
        if(expM2gtM1) {
            range *= -1;
        }
    }
    return(range);
}


inline int set_ISIZE(CMappingResult &m1, CMappingResult &m2, int readLength)
{
    bool expM2gtM1 = (m1.strand = '+');
    int sep = getSep(m1.uiPOS, m2.uiPOS, expM2gtM1);
    m2.ISIZE = m1.ISIZE = sep + readLength;
    return(sep);
}

bool CPairedReadsMapping::bookKeepValidMappedPairs
(CBestPairedMapping& bestMP, CPairedMappingFlags& flags, CMappingResult &m1, CMappingResult &m2)
{
    // Require reads mapped to the same ref sequence
    if (isValidPaired(m1, m2, opt)) {
        // TODO: double check if opt.readLength is good enough
        int sep =  set_ISIZE(m1, m2, (int)opt.readLength);
        if (opt.disLB <= sep && sep <= opt.disUB) {
            flags.bNoMapping = false;
            bestMP.update(m1, m2, opt.bExcludeAmbiguousReads);
            /*
            if (bPrintAllMapping) {
                this->printAMappedPair(table, m1, m2, );
            }*/
        } else {
            flags.sepMore |= (sep > opt.disUB);
            flags.sepLess |= (sep < opt.disLB);
        }
        bool firstEndFirst = m1.uiPOS > m2.uiPOS;
        flags.pairedOnExpStrand |= (isExpPairedMappedStrand(firstEndFirst, m1.strand, m2.strand, flags.mapSOLiDRead));
        return(true);
    }
    return(false);
}

int CPairedReadsMapping::printValidMappedPair(const GIndexTQ& table1, const GIndexTQ& table2, CMappingResult& m1, CMappingResult& m2, int validMappedPairNo)
{
    for (unsigned int i = 0; i < this->alignmentsQ[0].load; i++) {
        for (unsigned int j = 0; j < this->alignmentsQ[1].load; j++) {
            getSingleMappingIndex(*table1.pgenomeNT, this->alignmentsQ[0], i, m1);
            getSingleMappingIndex(*table2.pgenomeNT, this->alignmentsQ[1], j, m2);
            // Require reads mapped to the same ref sequence
            if (isValidPaired(m1, m2, opt)) {
                int sep =  set_ISIZE(m1, m2, (int)opt.readLength);
                if (opt.disLB <= sep && sep <= opt.disUB) {
                    this->printAMappedPair(table1, table2, m1, m2, validMappedPairNo);
                }
            }
        }
    }
    return(validMappedPairNo);
}

int CPairedReadsMapping::printBestMappedPair(const GIndexTQ& table1, const GIndexTQ& table2, CMappingResult& m1, CMappingResult& m2, int minMismatchNo, int bestMappedPairNo)
{
    for (unsigned int i = 0; i < this->alignmentsQ[0].load; i++) {
        for (unsigned int j = 0; j < this->alignmentsQ[1].load; j++) {
            getSingleMappingIndex(*table1.pgenomeNT, this->alignmentsQ[0], i, m1);
            getSingleMappingIndex(*table2.pgenomeNT, this->alignmentsQ[1], j, m2);
            // Require reads mapped to the same ref sequence
            if (isValidPaired(m1, m2, opt)) {
                int sep =  set_ISIZE(m1, m2, (int)opt.readLength);
                if (opt.disLB <= sep && sep <= opt.disUB) {
                    if ((int)(m1.uiDiff + m2.uiDiff) == minMismatchNo) {
                        this->printAMappedPair(table1, table2, m1, m2, bestMappedPairNo);
                    }
                }
            }
        }
    }
    return(bestMappedPairNo);
}

int CPairedReadsMapping::dealMappedPairedReads(const GIndexTQ& table1, const GIndexTQ& table2)
{
    CAlignmentsQ& q1 = this->alignmentsQ[0];
    CAlignmentsQ& q2 = this->alignmentsQ[1];
    CMappingResult m1(q1, opt.readLength);
    CMappingResult m2(q2, opt.readLength2);

    CBestPairedMapping bestMP;
    CPairedMappingFlags flags;
    flags.bPrintAllMapping = opt.bGetAllAlignments && !opt.bExcludeAmbiguousPaired;
    flags.mapSOLiDRead = this->opt.bMappedSOLiDRead; // Todo check if that selection is corrected

    bool samFormat = this->isSamFormat();
    if (!flags.mapSOLiDRead) {
        // get mapping info
        getReadQscores4Solexa(q1, m1, samFormat);
        getReadQscores4Solexa(q2, m2, samFormat);
    }

    for (unsigned int i = 0; i < q1.load; i++) {
        for (unsigned int j = 0; j < q2.load; j++) {
            getSingleMappingIndex(*table1.pgenomeNT, q1, i, m1);
            getSingleMappingIndex(*table2.pgenomeNT, q2, j, m2);
            bookKeepValidMappedPairs(bestMP, flags, m1, m2);
        }
    }
    this->dealMappingResult(table1, table2, m1, m2, bestMP, flags);
    return(0);
}

/*
 * This function do statistics and print mapped result.
 */
int CPairedReadsMapping::dealMappingResult(const GIndexTQ& table1, const GIndexTQ& table2, CMappingResult& m1, CMappingResult& m2, CBestPairedMapping& bestMP, CPairedMappingFlags& flags)
{
    bool samFormat = this->isSamFormat();
    if (flags.bNoMapping) {
        dealNoMapping(table1, table2, m1, m2);
        this->bookNoMappedKeepPairs(flags.sepMore, flags.sepLess, flags.pairedOnExpStrand);
        return(0);
    } else {
        this->getPairedRInfo(table1, table2, m1, m2, samFormat);
        if(flags.bPrintAllMapping) {
            this->printValidMappedPair(table1, table2, m1, m2, bestMP.validMappingNo);
        } else {
            this->dealBestMapping(table1, table2, bestMP, bestMP.bm1, bestMP.bm2);
        }
        this->bookKeepMappedPairs(bestMP);
        return(1);
    }
}

void CPairedReadsMapping::dealNoMapping(const GIndexTQ& table1, const GIndexTQ& table2, CMappingResult& m1, CMappingResult& m2)
{
    setSamFlags4OnlyOneEndMapped(m1, true);
    const int endId = 1;
    this->dealMappedSingleRead(table1, this->alignmentsQ[0], m1, endId == 1);
    setSamFlags4OnlyOneEndMapped(m2,false);
    this->dealMappedSingleRead(table2, this->alignmentsQ[1], m2, endId == 2);
}

int CPairedReadsMapping::dealBestMapping(const GIndexTQ& table1, const GIndexTQ& table2, CBestPairedMapping& bestMP, CMappingResult& m1, CMappingResult& m2)
{
    if(bestMP.validMappingNo == 1) { // -A -e & 1 mapping
        this->printAMappedPair(table1, table2, bestMP.bm1, bestMP.bm2, 1);
    } else if (opt.bExcludeAmbiguousPaired && !opt.bGetAllAlignments && bestMP.bestMappingNo == 1) { // -e 1 best mapping
        this->printAMappedPair(table1, table2, bestMP.bm1, bestMP.bm2, 1);
    } else if (!opt.bExcludeAmbiguousPaired && !opt.bGetAllAlignments) { // -B
        this->printBestMappedPair(table1, table2, m1, m2, bestMP.minDiff, bestMP.bestMappingNo);
    }
    return(bestMP.bestMappingNo);
}

inline void CPairedReadsMapping::getPairedRSeqAndQ(const GIndexTQ& table1, const GIndexTQ& table2, CMappingResult &m1, CMappingResult &m2, bool samFormat)
{
}

inline void CPairedReadsMapping::getPairedMappedRef4Solexa(const GIndexTQ& table1, const GIndexTQ& table2, CMappingResult &m1, CMappingResult &m2, bool samFormat)
{
    // get mapped reference sequence
    bool bNoRef = samFormat || !opt.bPrintRef4PairedInMapping;
    if(opt.bMappedLongRead) {
        getLongRefSeq(table1, m1, bNoRef);
    } else {
        getSingleMappingSeq4Solexa(table1, m1, bNoRef);
    }
    // TODO double check if opt.bMappedLongRead2 is set correctly.
    if(opt.bMappedLongRead2) {
        getLongRefSeq(table2, m2, bNoRef);
    } else {
        getSingleMappingSeq4Solexa(table2, m2, bNoRef);
    }

}

inline void CPairedReadsMapping::getMappedPairedRef(const GIndexTQ& table1, const GIndexTQ& table2, CMappingResult &m1, CMappingResult &m2, bool samFormat)
{
    // This should be called for every mapped paired
    if (this->opt.bPrintAlignments) {
        CAlignmentsQ& q1 = this->alignmentsQ[0];
        CAlignmentsQ& q2 = this->alignmentsQ[1];
        if (table1.bMapReadInColors) {
            getSingleMappingSeqAndQ4SOLiD(table1, q1, q1.read, m1, samFormat);
            getSingleMappingSeqAndQ4SOLiD(table2, q2, q1.read, m2, samFormat);
        } else {
            this->getPairedMappedRef4Solexa(table1, table2, m1, m2, samFormat);
        }
    }
}

// The function currnetly design to get the quality score for SOLEXA read only
inline void CPairedReadsMapping::getPairedRInfo(const GIndexTQ& table1, const GIndexTQ& table2, CMappingResult &m1, CMappingResult &m2, bool samFormat)
{
    if (this->opt.bPrintAlignments) {
        bool bMapReadInColors = table1.bMapReadInColors;
        if (bMapReadInColors) {
            // TODO: Distinguish if the following must be completed or it is not necessary 07/10.
            /*
            getSingleMappingSeqAndQ4SOLiD(table1, this->alignmentsQ[0], m1, samFormat);
            getSingleMappingSeqAndQ4SOLiD(table2, this->alignmentsQ[1], m2, samFormat);
            */
        } else {
            if(opt.bPrintPairedRQ) {
                // Get quality score
                getQscores4Solexa(this->alignmentsQ[0], m1, samFormat);
                getQscores4Solexa(this->alignmentsQ[1], m2, samFormat);
            }
            // this->getPairedMappedRef4Solexa(table, m1, m2, samFormat);
        }
    }
}

// F3read and R3read are index in AlignmentQ
void CPairedReadsMapping::printAMappedPair\
(const GIndexTQ& table1, const GIndexTQ& table2, CMappingResult &m1, CMappingResult &m2, int noPairedLoc)
{
    if (this->opt.bPrintAlignments) {
        set_ISIZE(m1, m2, (int)opt.readLength);
        bool samFormat = this->isSamFormat();
        // TODO get the mapped referenc sequence for every pair
        this->getMappedPairedRef(table1, table2, m1, m2, samFormat);
        if (samFormat) { // sam format
            printAPairedMappingInSam(this->AlignResult, m1, m2);
        } else {
            // string category = getCategory(m1.strand,m2.strand, m1.ISIZE, opt.disLB, opt.disUB);
            printAPairedMappingInPerM(this->AlignResult, m1, m2, noPairedLoc, opt.bPrintNM);
        }
    }
}

int CPairedReadsMapping::printLogFile(const char* inputFile)
{
    ofstream logFile(opt.logFileN, ofstream::app);
    const char* readSetName = getBasename(inputFile).c_str();

    if (logFile.good()) {
        this->printCommand(logFile, opt.fullCommand);
        this->printMappedPairStats(logFile, readSetName, opt.subDiffThreshold);
        this->printCommand(cout, opt.fullCommand);
        this->printMappedPairStats(cout, readSetName, opt.subDiffThreshold);
    } else {
        this->printCommand(cout, opt.fullCommand);
        this->printMappedPairStats(cout, readSetName, opt.subDiffThreshold);
    }
    logFile.close();
    return(0);
}

int CPairedReadsMapping::dealMappedSingleRead\
(const GIndexTQ& table, CAlignmentsQ &Q, CMappingResult &m, bool bFirstEnd)
{
    bool samFormat = this->isSamFormat();
    if (Q.load > 0 && samFormat) {
        for (unsigned int i = 0; i < Q.load && i < Q.iMaxCapacity; i++) {
            if (this->opt.bPrintAlignments) {
                getSingleMappingInfo(table, Q, i, m, samFormat, this->opt.bPrintNH);
                this->printSingleEndReads(m);
            }
        }
    } else if (this->opt.bPrintUnMappedReads) {
        // check if the seq is not decoded for SOLiD read?
        // check if code has never reach here?
        this->dealMissedRead(m);
    }
    if (Q.load > 0) {
        if (bFirstEnd) {
            this->noOfSingle1stEndMapped++;
        } else {
            this->noOfSingle2ndEndMapped++;
        }
    }
    return(this->iMissReadCounter++);
}

void CPairedReadsMapping::printMappedPairStats\
(ostream& out, const char* readSetName, unsigned int uiSubThreshold)
{
    // string readSetName = getBasename(readSet.readSetName());
    out << '\n';
    out << readSetName << ", #Pairs, " << this->iReadCounter << ", ";
    out << "#Mapped Pairs, " << this->noOfPairsInRange << ", ";
    out << "#Multi-mapped Pairs, " << this->noOfAmbiguousPairs << "\n";
    out << readSetName << ", ";
    out << "#Pairs sep more, " << this->noOfPairsSepMore << ", ";
    out << "#Pairs sep less, " << this->noOfPairsSepLess << ", ";
    out << "#Pairs sep more and less, " << this->noOfPairsSepMoreAndLess << endl;
    out << "#Pairs on exp strands " << this->noOfExpMappedPairedStrand << endl;
    out << readSetName << ", ";
    out << "#Mapped single 1st end, " << this->noOfSingle1stEndMapped << endl;
    out << "#Mapped single 2nd end, " << this->noOfSingle2ndEndMapped << endl;
    out << readSetName << ", ";
    unsigned int i;
    for (i = 0; i <= uiSubThreshold; i++) {
        out << "Sub" << i << ", " << iMapDiffCount[i] << ", " ;
    }
    out << endl;
}

/*
unsigned int CPairedReadsMapping::getPairedReadSetSize
(CReadInBitsSet& setA1, CReadInBitsSet& setA2, CReadInBitsSet& setB1, CReadInBitsSet& setB2)
{
    int setAsize = this->checkPairedReadSetSize(setA1, setA2);
    int setBsize = this->checkPairedReadSetSize(setB1, setB2);
    return(min(setAsize, setBsize));
}*/

int CPairedReadsMapping::dealMappedLongPairedRead
(const GIndexTQ& table1,  const GIndexTQ& table2,
 CMappingResult& m1, CMappingResult& m2)
{
    CAlignmentsQ& q1 = this->alignmentsQ[0];
    CAlignmentsQ& q2 = this->alignmentsQ[1];
    CBestPairedMapping bestMP;
    CPairedMappingFlags flags;
    flags.bPrintAllMapping = opt.bGetAllAlignments && !opt.bExcludeAmbiguousPaired;
    flags.mapSOLiDRead = false; // TODO: extend to SOLiD read

    for (unsigned int i = 0; i < q1.load; i++) {
        for (unsigned int j = 0; j < q2.load; j++) {
            // get mapping info
            getSingleMappingIndex(*table1.pgenomeNT, q1, i, m1);
            getSingleMappingIndex(*table2.pgenomeNT, q2, j, m2);
            bookKeepValidMappedPairs(bestMP, flags, m1, m2);
        }
    }
    return(this->dealMappingResult(table1, table2, m1, m2, bestMP, flags));
}

bool CPairedReadsMapping::isMappingSolidPaired(CLongReadsSet& set1, CReadInBitsSet& set2, const GIndexTQ& table1, const GIndexTQ& table2)
{
    if(opt.bMappedSOLiDRead && set1.isSOLiDReadFormat() && set2.isSOLiDReadFormat() \
            && table1.bMapReadInColors && table2.bMapReadInColors) {
        return(true);
    } else if(!opt.bMappedSOLiDRead && !set1.isSOLiDReadFormat() && !set2.isSOLiDReadFormat() \
              && !table1.bMapReadInColors && !table2.bMapReadInColors) {
        return(false);
    } else {
        LOG_INFO("\nInfo %d: Inconsistent index and input read (SOLiD vs Illumina)!", ERROR_LOG);
        return(false);
    }
}

